Material adapted from Grant R. McDermott’s notes
NOTE: the examples in this notebook focus only on vector-based spatial analysis. For information on raster-based spatial analysis, check some of the resources from Data Carpentry:
https://datacarpentry.org/organization-geospatial/
https://datacarpentry.org/r-raster-vector-geospatial/
(requirements vary by OS)
R provides bindings to powerful open-source GIS libraries. These include the - Geospatial Data Abstraction Library (GDAL) - Interface to Geometry Engine Open Source (GEOS) API suite - Access to projection and transformation operations from the PROJ.4 library.
sf, lwgeom, maps, mapdata, spData, tigris, tidycensus, tmap, tmaptoolsNote that you need ggplot2 version 3.0.0 or above for the required plotting support of sf objects.
library(tidyverse)
library(tidycensus)
library(sf)
library(tmap)
library(maps)
library(mapdata)
library(tigris)
We will access some data from the US Census Bureau through the tidycensus package. This will require a Census API key, which you can request here. Once that is done, you can set it using the tidycensus::census_api_key() function. You can use the "install = T" option to save your key.
Additional References: - Overview of Coordinate Reference Systems (CRS) in R
Spatial data (like all coordinate-based systems) only make sense relative to some fixed point. That fixed point is what the Coordinate Reference Systems, or CRS, is trying to set. In R, we define the CRS with a “proj4string”, which is based on the syntax of the PROJ.4 geospatial library.
In R, the notation used to describe the CRS is proj4string from the PROJ.4 library. It looks like this:
+init=epsg:4121 +proj=longlat +ellps=GRS80
+datum=GGRS87 +no_defs +towgs84=-199.87,74.79,246.62
Whenever we try to plot (some part of) the earth on a map, we are effectively trying to project a 3-D object onto a 2-D surface. This will necessarily create some kind of distortion. Different types of map projections limit distortions for some parts of the world at the expense of others. Distortions are particularly acute in the extreme latitudes for some projections, so try to choose according to the specific needs of your research question. For example, take a look at how badly the (in)famous Mercator projection distorts parts of the global map (source)
Animating the Mercator projection to the true size of each country in relation to all the others.
— Neil Kaye (@neilrkaye) October 12, 2018
Focusing on a single country helps to see effect best.#dataviz #maps #GIS #projectionmapping #mapping pic.twitter.com/clpCiluS1z
You can also “orient” your projection to a specific latitude and longitude using the PROJ.4 syntax.
sf packageR has long provided excellent support for spatial analysis and plotting (primarily through the sp, rgdal, rgeos, and raster packages). However, until recently, the complex structure of spatial data needed a set of equally complex spatial objects in R. The sf package provides simple features access for R. The “sf” stands for simple features, which is a simple standard for representing the spatial geometries of real-world objects on a computer. See the documentation and excellent sf vignettes for more details.
These objects — i.e. “features” — could include a tree, a building, a country’s border, or the entire globe. The point is that they are characterized by a common set of rules, defining everything from how they are stored on our computer to which geometrical operations can be applied them. Of greater importance to us here, however, is the fact that sf represents these features in R as data frames. This means that all of our favorite tidyverse operations can be applied to spatial data.
Most of the functions in the sf package start with the prefix st_. This stands for spatial and temporal, and a basic command of this package is easy enough once you remember that you are probably looking for st_SOMETHING().
As a quick demonstration we will read in the Florida counties shapefile. We use the st_read() command and the sf package will handle all the heavy lifting behind the scenes.
The shapefile was downloaded from: http://geodata.myflorida.com/datasets/swfwmd::florida-counties
# library(sf) ## Already loaded
## Location of our shapefile (here: bundled together with the sf package)
loc_counties <- "./Florida_Counties/Florida_Counties.shp"
## Read the shapefile into R
fla_counties <- st_read(loc_counties, quiet = TRUE)
Print out the fla object that we just created and take a look at its structure.
fla_counties
Simple feature collection with 67 features and 7 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -87.62601 ymin: 24.54522 xmax: -80.03095 ymax: 30.99702
CRS: 4326
First 10 features:
OBJECTID DEPCODE COUNTY COUNTYNAME DATESTAMP ShapeSTAre ShapeSTLen
1 1 21 041 GILCHRIST 2000-05-16T00:00:00.000Z 9908353355 487300.0
2 2 54 107 PUTNAM 2000-05-16T00:00:00.000Z 23058687376 762967.7
3 3 62 123 TAYLOR 2000-05-16T00:00:00.000Z 28917470565 877252.7
4 4 46 091 OKALOOSA 2000-05-16T00:00:00.000Z 25621585036 1087058.0
5 5 7 013 CALHOUN 2000-05-16T00:00:00.000Z 16048085397 631344.0
6 6 50 099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969 1522726.9
7 7 56 111 ST. LUCIE 2000-05-16T00:00:00.000Z 16149484821 829518.7
8 8 51 101 PASCO 2000-05-16T00:00:00.000Z 21260807286 739390.6
9 9 20 039 GADSDEN 2000-05-16T00:00:00.000Z 14744963472 673612.4
10 10 37 073 LEON 2000-05-16T00:00:00.000Z 19584016138 807989.7
geometry
1 MULTIPOLYGON (((-82.65814 2...
2 MULTIPOLYGON (((-81.58084 2...
3 MULTIPOLYGON (((-83.73037 3...
4 MULTIPOLYGON (((-86.39159 3...
5 MULTIPOLYGON (((-84.93266 3...
6 MULTIPOLYGON (((-80.1119 26...
7 MULTIPOLYGON (((-80.23992 2...
8 MULTIPOLYGON (((-82.05534 2...
9 MULTIPOLYGON (((-84.28258 3...
10 MULTIPOLYGON (((-84.00743 3...
The object has the familar tibble-style output that we are used to (e.g. it only prints the first 10 rows of the data). However, it also has some additional information in the header, like a description of the geometry type (MULTIPOLYGON). Notice the geometry column right at the end of the data frame. This geometry column is how the sf package achieves much of its magic: It stores the geometries of each row element in its own list column. For example, we could print out the coordinates needed to plot the first element in our data frame, by typing fla_counties$geometry[[1]]. In contrast, compare with how complicated the structure of a traditional spatial object is by running, say, str(as(fla, "Spatial")).]
# this may print long output
# str(as(fla, "Spatial"))
ggplot2Plotting sf objects is incredibly easy thanks to the package’s integration with ggplot2. Plotting sf maps using the base R plot function is also easy (and often faster). The key geom to remember is geom_sf(). For example:
# library(tidyverse) ## Already loaded
ggplot(data = fla_counties) +
geom_sf(aes(fill = ShapeSTAre), alpha=0.8, col="white") +
scale_fill_viridis_c(name = "Area") +
ggtitle("Counties of Florida")
Largest counties in Florida:
| Rank | Land Area | County |
|---|---|---|
| 1 | 1,998.32 sq mi | Collier, FL |
| 2 | 1,969.76 sq mi | Palm Beach, FL |
| 3 | 1,897.71 sq mi | Miami-Dade, FL |
| 4 | 1,797.84 sq mi | Polk, FL |
fla_counties %>%
top_n(ShapeSTAre, n = 4)
Simple feature collection with 4 features and 7 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -82.10574 ymin: 25.14405 xmax: -80.03095 ymax: 28.36218
CRS: 4326
OBJECTID DEPCODE COUNTY COUNTYNAME DATESTAMP ShapeSTAre ShapeSTLen
1 6 50 099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969 1522727
2 43 53 105 POLK 2000-05-16T00:00:00.000Z 56039729500 1273588
3 47 11 021 COLLIER 2000-05-16T00:00:00.000Z 55941206046 1903355
4 48 13 025 DADE 2000-05-16T00:00:00.000Z 55073375764 1774363
geometry
1 MULTIPOLYGON (((-80.1119 26...
2 MULTIPOLYGON (((-81.65737 2...
3 MULTIPOLYGON (((-81.27162 2...
4 MULTIPOLYGON (((-80.87314 2...
Florida Lakes
Obtained from https://geodata.dep.state.fl.us/datasets/97b765ff2b70400d8bcab23fbe2a5e88_0
loc_lakes <- "./Florida_Lakes/Florida_Lakes.shp"
## Read the shapefile into R
fla_lakes <- st_read(loc_lakes, quiet = TRUE)
Plot lakes in Polk county
dim(fla_lakes)
[1] 4243 7
colnames(fla_lakes)
[1] "PERIMETER" "NAME" "COUNTY" "OBJECTID" "SHAPEAREA" "SHAPELEN" "geometry"
fla_lakes %>%
filter(COUNTY == "POLK") %>%
ggplot() +
geom_sf(aes(fill = SHAPEAREA), alpha=0.8, col="white") +
#scale_fill_viridis_c(name = "Area") +
ggtitle("Lakes in Polk County")
Biggest lakes in Polk County
fla_lakes %>%
filter(COUNTY == "POLK") %>%
top_n(SHAPEAREA, n = 10)
Simple feature collection with 10 features and 6 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -81.951 ymin: 27.66273 xmax: -81.36433 ymax: 28.10574
CRS: 4326
PERIMETER NAME COUNTY OBJECTID SHAPEAREA SHAPELEN geometry
1 16855.88 Lake Hamilton POLK 67 8668138 16855.88 MULTIPOLYGON (((-81.65413 2...
2 16652.12 Reedy Lake POLK 76 14191962 16652.12 MULTIPOLYGON (((-81.52262 2...
3 21815.34 Lake Parker POLK 1608 8533253 21815.34 MULTIPOLYGON (((-81.93106 2...
4 21007.95 Lake Marion POLK 1936 12226186 21007.95 MULTIPOLYGON (((-81.53609 2...
5 31582.90 Lake Pierce POLK 1937 15604508 31582.90 MULTIPOLYGON (((-81.49978 2...
6 34664.94 Crooked Lake POLK 2690 16885793 34664.94 MULTIPOLYGON (((-81.57241 2...
7 21529.93 Lake Weohyakapka POLK 2700 30549713 21529.93 MULTIPOLYGON (((-81.4103 27...
8 18557.39 Lake Rosalie POLK 2701 18358611 18557.39 MULTIPOLYGON (((-81.39447 2...
9 19950.41 Lake Hancock POLK 2703 18549186 19950.40 MULTIPOLYGON (((-81.83617 2...
10 20955.50 Lake Arbuckle POLK 3907 15412988 20955.50 MULTIPOLYGON (((-81.41045 2...
To reproject an sf object to a different CRS, we can use sf::st_transform().
fla_counties %>%
st_transform(crs = "+proj=moll +ellps=WGS84") %>% ## Reprojecting to a Mollweide CRS
head(2) ## Saving vertical space
Simple feature collection with 2 features and 7 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -7625474 ymin: 3564689 xmax: -7485526 ymax: 3636317
CRS: +proj=moll +ellps=WGS84
OBJECTID DEPCODE COUNTY COUNTYNAME DATESTAMP ShapeSTAre ShapeSTLen
1 1 21 041 GILCHRIST 2000-05-16T00:00:00.000Z 9908353355 487300.0
2 2 54 107 PUTNAM 2000-05-16T00:00:00.000Z 23058687376 762967.7
geometry
1 MULTIPOLYGON (((-7586204 36...
2 MULTIPOLYGON (((-7486900 36...
Or, we can specify a common projection directly in the ggplot() call using coord_sf(). This is often the most convenient approach when you are combining multiple sf data frames in the same plot.
ggplot(data = fla_counties) +
geom_sf(aes(fill = ShapeSTAre), alpha=0.8, col="white") +
scale_fill_viridis_c(name = "Area") +
coord_sf(crs = "+proj=moll +ellps=WGS84") +
labs(
title = "Counties of Florida",
subtitle = "Mollweide projection"
)
dplyr and tidyrThe tidyverse approach to data wrangling carries over very smoothly to sf objects. For example, the standard dplyr verbs like filter(), mutate() and select() all work. You can also perform group_by() and summarise() operations as usual. Furthermore, the dplyr family of join functions also work, which can be especially handy when combining different datasets by (say) Federal Information Processing Standard (FIPS) code or some other attribute. However, this presumes that only one of the objects has a specialized geometry column. In other words, it works when you are joining an sf object with a normal data frame. In cases where you want to join two sf objects based on their geometries, there is a specialized st_join() function.
Alongside all the tidyverse functionality, the sf package comes with a full suite of geometrical operations. You should take a look at the third sf vignette or the Geocomputation with R book to get a complete overview.
So-called unary operations are applied to a single object. For instance, you can “melt” sub-elements of an sf object (e.g. counties) into larger elements (e.g. states) using sf::st_union():
fla_counties %>%
st_union() %>%
ggplot() +
geom_sf(fill=NA, col="black") +
labs(title = "Outline of Florida")
Or, you can get the st_area(), st_centroid(), st_boundary(), st_buffer(), etc. of an object using the appropriate command. For example:
fla_counties %>% st_area() %>% head(5) ## Only show the area of the first five counties to save space.
Units: [m^2]
[1] 920490582 2142439339 2685252824 2369242717 1487627072
And:
fla_centroid <- st_centroid(fla_counties)
ggplot(fla_counties) +
geom_sf(fill = "black", alpha = 0.8, col = "white") +
geom_sf(data = fla_centroid, col = "red") + ## Notice how easy it is to combine different sf objects
labs(
title = "Counties of Florida",
subtitle = "Centroids in red"
)
Another set of so-called binary operations can be applied to multiple objects. So, we can get things like the distance between two spatial objects using sf::st_distance(). In the below example, we get the distance from Polk county to Palm Beach county:
fla_counties %>%
filter(COUNTYNAME %in% c("POLK", "PALM BEACH"))
Simple feature collection with 2 features and 7 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -82.10574 ymin: 26.32093 xmax: -80.03095 ymax: 28.36218
CRS: 4326
OBJECTID DEPCODE COUNTY COUNTYNAME DATESTAMP ShapeSTAre ShapeSTLen
1 6 50 099 PALM BEACH 2000-05-16T00:00:00.000Z 61491591969 1522727
2 43 53 105 POLK 2000-05-16T00:00:00.000Z 56039729500 1273588
geometry
1 MULTIPOLYGON (((-80.1119 26...
2 MULTIPOLYGON (((-81.65737 2...
polk_palmbeach <- fla_counties %>%
filter(COUNTYNAME %in% c("POLK", "PALM BEACH"))
palmbeach <- fla_counties %>%
filter(COUNTYNAME %in% c("PALM BEACH"))
## Use "by_element = T" to give a vector instead of the default pairwise matrix
ab_dist <- st_distance(polk_palmbeach, palmbeach, by_element = T)
## We can use the `units` package (already installed as sf dependency) to convert to kilometres
ab_dist <- ab_dist %>% units::set_units(km) %>% round()
Create plot:
ggplot(fla_counties) +
geom_sf(fill = "black", alpha = 0.8, col = "white") +
geom_sf(data = fla_counties %>% filter(COUNTYNAME %in% c("POLK", "PALM BEACH")),
aes(fill = COUNTYNAME), col = "white") +
labs(
title = "Calculating distances",
subtitle = paste0("The distance between Polk and Palm Beach is ", ab_dist[2], " km")
) +
theme(legend.title = element_blank())
We can easily import external shapefiles, KML files, etc., into R using the generic sf::st_read() function.
However, R also provides access to a large number of base maps — e.g. countries of the world, US states and counties, etc. — through the maps, (higher resolution) mapdata and spData packages, as well as a whole ecosystem of more specialized GIS libraries.
# library(maps) ## Already loaded
world1 <- st_as_sf(map("world", plot = FALSE, fill = TRUE))
world1 %>%
ggplot() +
geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
labs(title = "The world", subtitle = "Mercator projection")
All of the usual sf functions and transformations can then be applied. For example, we can reproject the above world map onto the Lambert Azimuthal Equal Area projection (and further orientate it at the South Pole) as follows.
world2 <-
st_transform(
world1,
"+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs"
)
world2 %>%
ggplot() +
geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
labs(title = "The world", subtitle = "Lambert Azimuthal Equal Area projection")
As we have already seen, most map projections work great “out of the box” with sf. One notable exception is the Winkel tripel projection. This is the preferred global map projection of National Geographic and requires a bit more work to get it to play nicely with sf and ggplot2 (as detailed in this thread). Here’s a quick example of how to do it:
# library(lwgeom) ## Already loaded
world3 <- lwgeom::st_transform_proj(world1, crs = "+proj=wintri +datum=WGS84 +no_defs +over")
## Don't necessarily need a graticule, but if you do then define it manually:
gr <-
st_graticule(lat = c(-89.9,seq(-80,80,20),89.9)) %>%
lwgeom::st_transform_proj(crs = "+proj=wintri +datum=WGS84 +no_defs +over")
world3 %>%
ggplot() +
geom_sf(data = gr, color = "#cccccc", size = 0.15) + ## Manual graticule
geom_sf(fill = "grey80", col = "grey40", lwd = 0.3) +
coord_sf(datum = NA) +
labs(title = "The world", subtitle = "Winkel tripel projection")
The maps and mapdata packages have detailed county- and province-level data for several individual nations. We can still use it to extract a specific country’s border using some intuitive syntax. For example, we could plot a base map of Colombia as follows.
colombia <- st_as_sf(map("world", "colombia", plot = FALSE, fill = TRUE))
colombia %>%
ggplot() +
geom_sf(fill="blue", col=NA)
Aside: You can detach the maps package once you are finished using it, since it avoids potential namespace conflicts with purrr::map.
detach(package:maps) ## To avoid potential purrr::map() conflicts
tigris and tidycensusWorking with Census data may require some extra work. You need to register on the website, then download data from various years or geographies separately, merge these individual files, etc. Thankfully, this has recently become much easier thanks to the Census API and — for R at least — the tigris and tidycensus packages from Kyle Walker. Check his tutorial on his website.
NOTE: Before continuing with this section, you will first need to request an API key from the Census. Once that’s done, you can set it using the
tidycensus::census_api_key()function. I recommend using the “install = T” option to save your key for future use.
We start by loading the packages and setting some options for optimized use with the sf package.
# library(tidycensus) ## Already loaded
# library(tigris) ## Already loaded
options(tigris_class = "sf")
options(tigris_use_cache = TRUE)
# census_api_key("YOUR KEY HERE", install = T)
Let’s say that our goal is to provide a snapshot of Census rental estimates across different cities in Florida. We start by downloading tract-level rental data for Florida using the tidycensus::get_acs() function. Note that you will need to look up the correct ID variable (in this case: “DP04_0134”).
rent <-
tidycensus::get_acs(
geography = "tract", variables = "DP04_0134",
state = c("FL"), geometry = TRUE
)
Getting data from the 2014-2018 5-year ACS
Using the ACS Data Profile
Using FIPS code '12' for state 'FL'
rent
Simple feature collection with 4245 features and 5 fields (with 33 geometries empty)
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -87.63494 ymin: 24.5231 xmax: -80.03136 ymax: 31.00089
CRS: 4269
First 10 features:
GEOID NAME variable estimate moe
1 12001000200 Census Tract 2, Alachua County, Florida DP04_0134 853 74
2 12001000301 Census Tract 3.01, Alachua County, Florida DP04_0134 840 105
3 12001000302 Census Tract 3.02, Alachua County, Florida DP04_0134 933 92
4 12001000400 Census Tract 4, Alachua County, Florida DP04_0134 1019 432
5 12001000500 Census Tract 5, Alachua County, Florida DP04_0134 861 90
6 12001000600 Census Tract 6, Alachua County, Florida DP04_0134 865 169
7 12001000700 Census Tract 7, Alachua County, Florida DP04_0134 835 41
8 12001000806 Census Tract 8.06, Alachua County, Florida DP04_0134 786 54
9 12001000808 Census Tract 8.08, Alachua County, Florida DP04_0134 906 43
10 12001000809 Census Tract 8.09, Alachua County, Florida DP04_0134 893 48
geometry
1 MULTIPOLYGON (((-82.33935 2...
2 MULTIPOLYGON (((-82.33922 2...
3 MULTIPOLYGON (((-82.33916 2...
4 MULTIPOLYGON (((-82.32413 2...
5 MULTIPOLYGON (((-82.33082 2...
6 MULTIPOLYGON (((-82.31147 2...
7 MULTIPOLYGON (((-82.32629 2...
8 MULTIPOLYGON (((-82.33938 2...
9 MULTIPOLYGON (((-82.34001 2...
10 MULTIPOLYGON (((-82.35586 2...
This returns an sf object, which we can plot directly.
rent %>%
ggplot() +
geom_sf(aes(fill = estimate, color = estimate)) +
scale_fill_viridis_c(name = "Rent ($)", labels = scales::comma) +
scale_color_viridis_c(name = "Rent ($)", labels = scales::comma) +
labs(
title = "Rental rates across Florida",
caption = "Data: US Census Bureau"
)
The above map is very detailed. Perhaps we are not interested in the entire set of tract level data, but would rather get a sense of rents within some well-defined metropolitan areas? The tigris package has you covered here. For example, let’s say we want to compare average rents across three metros: Orlando, Tampa and Miami.
fl_metros <-
tigris::core_based_statistical_areas(cb = TRUE) %>%
# filter(GEOID %in% c("53000", "71000", "45000")) %>% ## Could use GEOIDs directly if you know them
filter(grepl("Orlando|Tampa|Miami", NAME)) %>%
select(metro_name = NAME)
Now we do a spatial join on our two data sets using the sf::st_join() function.
fl_rent <-
st_join(
rent,
fl_metros,
join = st_within, left = FALSE
)
fl_rent
Simple feature collection with 2349 features and 6 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -82.85243 ymin: 25.13807 xmax: -80.03136 ymax: 29.27677
CRS: 4269
First 10 features:
GEOID NAME variable estimate moe
223 12011010102 Census Tract 101.02, Broward County, Florida DP04_0134 1728 393
224 12011010103 Census Tract 101.03, Broward County, Florida DP04_0134 1138 138
225 12011010104 Census Tract 101.04, Broward County, Florida DP04_0134 1341 180
226 12011010200 Census Tract 102, Broward County, Florida DP04_0134 1262 79
227 12011010304 Census Tract 103.04, Broward County, Florida DP04_0134 1074 199
228 12011010305 Census Tract 103.05, Broward County, Florida DP04_0134 1146 45
229 12011010306 Census Tract 103.06, Broward County, Florida DP04_0134 1042 98
230 12011010307 Census Tract 103.07, Broward County, Florida DP04_0134 1142 59
231 12011010308 Census Tract 103.08, Broward County, Florida DP04_0134 1940 119
232 12011010401 Census Tract 104.01, Broward County, Florida DP04_0134 1732 102
metro_name geometry
223 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.09418 2...
224 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.09123 2...
225 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.08109 2...
226 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.10748 2...
227 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.11944 2...
228 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.12199 2...
229 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.10789 2...
230 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.12833 2...
231 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.15345 2...
232 Miami-Fort Lauderdale-West Palm Beach, FL MULTIPOLYGON (((-80.17015 2...
One useful way to summarize this data and compare across metros is with a histogram. Note that “regular” ggplot2 geoms and functions play perfectly nicely with sf objects (i.e. we are not limited to geom_sf()).
fl_rent %>%
ggplot(aes(x = estimate)) +
geom_histogram() +
facet_wrap(~metro_name, labeller = label_wrap_gen(width=10)) +
theme_bw()